function [profit] = calc_pi(vec,params)
%calculate profit in plants in space model
%   params has epsilon, theta, xi, q, sigma, grid_sz, plant_locs, D, BoW, R
% unpack parameters first, unless it's 0 locations open
% vec should be vector of decisions for plants
N_init=sum(vec);

if N_init==0
    profit=0;
else
    epsilon=params.epsilon;
    theta=params.theta;
    constant=(1/(epsilon-1))*((epsilon/(epsilon-1))^(-epsilon))*gamma(1-((epsilon-1)/theta));
    xi=params.xi;
    q=params.q; % if too high, profits go to infinity
    sigma=params.sigma; % if Zprod uses form Z=q/N^sigma
    delta=params.delta;
    % sigma=.060368; %if Zprod uses form q*exp(-N*sigma)

    grid_sz=params.grid_sz; % # of points to make in square, grid_sz-1=dimension
    Ds=params.D; %scalar, demand, could make a matrix for each s?

%     plants=zeros(N_init,1);
    BoW=zeros(N_init,1);
    R=zeros(N_init,1);
    k=1;
    for l=1:length(vec)
        if vec(l)==1
            plant_locs(k,:)=params.plant_locs(l,:);
            BoW(k)=params.BoW(l);
            R(k)=params.R(l);
            k=k+1;
        end
    end
    profit=0;
    for i=0:grid_sz
        for j=0:grid_sz
            s=[i/grid_sz j/grid_sz];
            D=Ds(i+1, j+1);
            inner_sum=0;
            for n=1:N_init
                plant=plant_locs(n,:);
                inner_sum=inner_sum+((BoW(n)*Zprod(q,N_init*delta^2,sigma)/(((norm(s-plant))/delta)+1))^theta);
            end
            profit=profit+D*inner_sum^((epsilon-1)/theta);
        end
    end
    profit=constant*profit-sum(R)*xi;
    profit=round(profit,10);
end

end

